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ABSTRACT 



A new technique for constructing "computational molecules" 
for linear finite difference operators is developed. The basic 
approach is one of approximating a two dimensional surface with 
a geometrically consistent interpolating polynomial of degree 
four or five. The desired finite difference operator is then 
developed from the polynomial. The resulting molecules are 
geometrically consistent and may be used to solve boundary 
value problems without the use of fictitious points. 

Molecules for the biharmonic operator with various 
boundary conditions included are presented in this paper, as 
well as molecules representing the boundary conditions for 
shear and moment along the free edge of a plate. 

The integrity of the molecules presented is proven by 
comparison of solutions for flat plate bending problems by 
finite difference with exact solutions from the literature. 
Convergence plots for each problem are also presented. 
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NOMENCLATURE 



A , ab 


Differentiation of some function A with respect 
to variables (a) and (b) . 


C i 


Constants for an interpolating polynomial of degree 
four. 


D 


p-t 3 

Flexural rigidity: £,T 

12 (1-y 2 ) 


D. 

l 


Constants for an interpolating polynomial of degree 


E 


Young's modulus of elasticity. 


E t 


Error of approximation. 


(x,y) 


Expression for variable terms of an interpolating 
polynomial of degree four. 


6<i) (x,y) 


Expression for variable terms of an interpolating 
polynomial of degree five. 


h 


Mesh size spacing. 


L 


Length or width. 


Mn 


Moment in the normal direction along the edge of a 
flat plate. 


N 


Number of divisions on one side of a square finite 
difference mesh. 


n 


Normal co-ordinate . 


P 


Applied pressure 


P(x,y) 


Interpolating polynomial. 


P i 


Values of polynomial at various nodal points. 


t 


Tangential co-ordinate. 


V 


Transverse shear along a free edge of a plate. 


X 


Horizontal co-ordinate. 


y 


Vertical co-ordinate 


6 


Plate deflection. 


u 


Poisson's ratio 



5 



{ } Matrix notation for column vector. 

( > Matrix notation for row vector. 

[ ] Square or rectangular matrix. 

Harmonic operator. 

Biharmonic operator. 

T Thickness of a thin plate. 

O' Order of 
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I . INTRODUCTION 



A. GENERAL DESCRIPTION 

Solving partial differential equations and boundary value 
problems by approximate numerical methods is not new. Such 
numerical solutions have become very popular in the past twenty 
years due to the development of the high speed digital computer. 

One widely used numerical method for solving partial 
differential equations is that of finite differences. The 
solution of equations by this method leads to the development 
of "computational molecules". Various techniques have been 
used to develope such molecules and numerous molecules are 
available in the literature (Refs. 1,2,3, &4) . 

Often these molecules tend to be geometrically inconsistent. 
This inconsistency results from combining molecules for the 
elementary differential operators, developed from different 
interpolating polynomials, in order to obtain a more complex 
operator. For instance, the molecule for the biharmonic operator 
as given by Salvadori and Baron (Ref. 1) utilizes two different 
interpolating formulas. The fourth derivatives are obtained 
from a polynomial of degree four. The mixed partial derivative 
is obtained from a bi-quadratic interpolating surface. 

Furthermore, the application of these molecules in the 
solution of boundary value problems requires the use of ficti- 
tious nodal points outside the physical boundaries of the 
problem. Relationships between these fictitious nodal points 
and those nodal points representing the physical problem must 
be pre-determined . When the fictitious points cannot be 

determined uniquely, the method fails. 
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Considering the above disadvantages of the classical finite 
difference operators found in the literature, the objective 
of the author's research was to develope a new technique for 
constructing linear finite difference operators. This new 
technique utilizes one consistent form of an interpolating 
polynomial to approximate a two dimensional surface from which 
molecules for all operators are constructed. Furthermore, 
basic operators incorporating physical boundary conditions 
are constructed using this technique. The use of such operators 
alleviates the need for fictitious nodal points. 

The technique proposed in this research is applicable 
for the construction of any finite difference operator. How- 
ever, the field of investigation for the purpose of this 
research was limited to the construction of those operators 
required for the solution of plate bending problems with various 
boundary conditions. The solutions of such problems will be 
used to justify the integrity of the constructed finite dif- 
ference operators . 

B. BASIC TECHNIQUE 

The basic technique employed to construct the molecules 
for a prescribed linear finite difference operator is one of 
approximating a two dimensional surface by means of a biquad- 
ratic or biquintic interpolating polynomial. The constants 
for such an interpolating polynomial are determined for each 
desired surface, which may include specified boundary condi- 
tions, using a finite difference grid. The polynomial is 
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then operated on with the desired partial differential 
operator to be approximated. 

The final result is a numerical approximation of the 
operator in terms of a computational molecule, grid size, 
nodal deflections, and nodal boundary condition values. 

All calculations required for the construction of these 
molecules were performed in double precision using an 
IBM-360 computer. 
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II. DERIVATION OF DIFFERENTIAL OPERATORS 



A. INTERPOLATING POLYNOMIALS 



The interpolating polynomials used to approximate a two 
dimensional surface for the development of a finite difference 
biharmonic operator in this work were either 



P (x,y) = C 1 y 4 x 4 + C 2 y 3 x 4 + C 3 y 2 x 4 + C 4 yx 4 + C 5 x 4 

+ Cgy 4 x 3 + C 7 y 3 x 3 + Cgy 2 x 3 + C^yx-^+C^gX 3 
+C 11 y4x2+C 12 y3 x 2+C 13 y2x2+C ]L4 yx2+c 15 x2 

+C 16 y4x +C 17 y ^ x +C 18 y2x +C 19 yx +C 20 X 
+C21Y 4 +c 22Y 3 +c 23Y 2 +c 24Y +c 25 



and in matrix notation 



(i) 



P (x,y) = < F (x,y) ) {C.:} • , 2 

JL — -L f f • • • f 



25 



or 



(2.1a) 



P (x,y) 



D-^y 3 x 3 + D 2 y 4 x 5 + D 3 y 3 x 3 + D 4 y 2 x 3 + D^yx^+ D^x 3 
+ D^y^x 4 + Dgy^x 4 + D^y ^x 4 +D^gy ^x 4 +D^^yx 4 +D^ 2 x 4 

+D 13 y ^ x2+D 14 y4x2+D 15 y2x2+D 16 y2x2+D 17 yx2+D 18 x2 

+D 19 y ^ x2+D 20 y4x2+D 21 y2x2+D 22 y2x2+D 23 yx2+D 24 x2 
+D 25 y5 x +D 26 y 4 x +D 2? y 3 x +D 2g y 2 x +D 29 yx +D 3Q x 

+d 3 iY 5 +d 32 y 4 +d 33 y 3 +d 34 y 2 +d 35 y +d 36 



( 2 . 2 ) 
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and in matrix notation 



P (x,y) = <G^(x,y)> {D. } 

1 i=l,2,...,36 (2.2a) 

Polynomials (2.1) (2.2) will remain polynomials of degree 

four or five respectively for both constant values of x and 
y. All partial derivatives with respect to x and y of order 
four exist. Either polynomial will consistently approximate 
a two dimensional surface. 

In certain cases, a fifth order polynomial such as (2.2) 
which will approximate a thirty-six nodal point grid surface 
is required in the construction of the biharmonic operator 
molecule with included boundary conditions. This requirement 
results from the need for keeping the error of approximation 
for a given molecule of an order equal to or greater than the 
square of the grid size spacing. The subject of approximation 
error will be further discussed later. 

B „ BASIC FINITE DIFFERENCE BIHARMONIC OPERATOR 

For the solution of a problem governed by a partial 
differential equation, in which the biharmonic operator appears, 
a numerical approximation of the biharmonic operator without 
boundary conditions must be constructed first. To construct 
the molecule for such an operator "'olynomial (2.1) was used to 
approximate a two dimensional s' ,:e defined by the ordinates 

at the nodes of a grid shown in ,ure (1) . 

To determine the unknown constants for the interpolating 
polynomial , the polynomial was evaluated at each of the 
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twenty-five nodal points In terms of the co-ordinate values 
of x and y. In matrix notation, the resulting system of 
equations are expressed as 



{p ± } 


= [M] {C-}. 

X 1 


where 




{p i } 


= (P(x i ,y i )} 


[M] 


1! 

— 1 

/s 

p* 



25 



(2.3) 



25 



< F (l) (x 2 / y 2 ) > 



< fU) ( x 25^25 } > 

i=l , 2 , . . . ,25_ 

Solving equation (2.3) for the unknown coefficients, yields 

(Ci) = [M] {P i } i=l,2, . . . ,25 (2.4) 

with which the interpolating polynomial is rewritten as 

P(x,y) = <F (i) (x,y)> CM]” 1 {Pi) iMl(2 25 (2 . 5) 

Operating on equation (2.5) with the biharmonic operator 
results in the following equation 

V 4 P(x,y) = < V 4 F (i) (x,y) > [M]" 1 {P ± }. 

x 1=1 , 2 , . . . , 25 (2.6) 

The quantity ( V 4 F^(x,y)) is now evaluated for a particular 
nodal point say m (x=2h,y=2h) . The resulting quantity 
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( (2h,2h)> [M] after numerical evaluation yields a 

five by five computational molecule for the finite difference 
biharmonic operator at nodal point (m) . 

The above molecule allows for the numerical approximation 
of the biharmonic operator in terms of known constants and 
nodal deflections. Evaluation of a computational molecule 
for a nodal point other than the center one is performed in 
a similar manner except that polynomial (2.2) and a six by 
six grid is used in order to achieve an acceptable bound on 
the error of approximation. 

C. FINITE DIFFERENCE BIHARMONIC OPERATORS WITH BOUNDARY 

CONDITIONS SPECIFIED 

To complete the solution of a boundary value problem, the 
applicable boundary conditions must be considered. If geo- 
metric considerations alone are required for the boundary 
conditions, a general finite difference molecule for the 
biharmonic operator incorporating the desired boundary con- 
ditions may be developed. The technique uses polynomial (2.2) 
to approximate a surface having the desired boundary conditions. 

1 . Operators With First Partial Derivative Boundary 
Conditions 

For the construction of a finite difference bihar- 
monic operator with a first partial with respect to y boundary 
condition, polynomial (2.2) is used to approximate a surface 
such as the one represented by Figure (2) where the nodal 
points along the top edge are repeated to accommodate both 
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values of deflection and the specified boundary conditions 
along that edge. 

The unknown coefficients for the interpolating 
polynomial are determined by evaluating the polynomial at each 
of the thirty distinct nodal points and the first partial with 
respect to y of the polynomial at each of the nodal points 
along the top edge. The resulting equations are 



P (x,y) = < G (x,y) > 


{D i } i=l,2, ... ,36 


(2.7) 


P v (x,y) = < G ^ (x,y) > 

f Y 


{D i } i=l,2, . . . ,36 


(2.8) 



The complete system of equations becomes 



< p i> 



(Mil (Di) 

X -i- / ^ / • * * / 



36 



(2.9) 



where 



{P.} 

1 



< P(xj ,yj ) j P, y (x k ,y k ) > 



[M^] 



< gU S.V> j=l,2 



30 



< G , ( y ) <x k ,y k )> 



k=31 , 32 , ... ,36 



i=l , 2 , . . . , 36 



and solving for the unknown coefficients 



{ Di } 



[M x ] 1 (Pi) 



r ”■ 1,2, . . . ,36 



( 2 . 10 ) 
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from which the interpolating polynomial may be written as 



P(x,y) = (G^Ufy)) [M -.]” 1 (P - } 

i=l , 2 , . . . ,36 (2.11) 

Operating on equation (2.11) with the biharmonic operator 
results in the following equation 



V 4 P(x,y) = < V 4 G (l) (x,y) > [M^ -1 (P^ 



i=l,2,...,36 (2.12) 



The quantity (V 4 G^(x,y)} is evaluated for any 
desired nodal point within the grid and the resulting numerical 
evaluation of (V 4 G^ (x,y) ) [M-jJ gives a six by six molecule 
for the biharmonic operator with the prescribed boundary 
conditions . 



For first partial derivatives with respect to x and 
to y along two perpendicular edges, the respective biharmonic 
molecule is developed using polynomial (2.2) to approximate 
the surface defined by ordinates at the nodes of the grid 
given in Figure (3) . The nodal points along two perpendicular 
edges are repeated to account for both the values of deflections 
and the boundary conditions at these nodal points. 

The coefficients for the interpolating polynomial are 
determined from the following relationships 



P(x,y) 

P y (x,y) 



(i) 



(i) 

,y 



(x,y) > 


{ d l } 


i=l , 2 , ... ,36 


(x,y) > 


{ Di } 




i=l , 2 , . . . , 36 



(2.13) 

(2.14) 
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(2.15) 



P X <*/Y> 



< (x,y) > {D^} 

f i=l , 2 , . , . , 36 



P (x,y) = 
,xy 



(i) 



G ;^(x,y)> 



( Di ) 



i=l , 2 , ,36 



(2.16) 



The mixed second partial was used in the development in order 
to have sufficient conditions to evaluate the thirty six 
coefficients of the interpolating polynomial and to provide 
for continuity at the common corner nodal point. 

The complete system of equations becomes 



{P,} = [M 2 ] {d.} 

i=l , 2 ,.. . ,36 



(2.17) 



where 



{pp = < p < x 1 'yj)! p / y< x k'yk>! p , x < x t'yp| p ,xy (x m'ym ) > 






G (l) (Xj ,yj) > 

J -L r £ f * » o f 



25 



< G / ( y ) (x k' y k } > k= 



k=26 , 27 , ... ,30 



< G( ^ (x il' Y £ } > 

' x * * £=31,32, . . . ,35 



^ G ,xy ^ x m' y m^ ^ m =36 



i=l , 2 , . . . , 36 
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Solving for the constants, {Dj_}, and operating on 
the resulting polynomial with the biharmonic operator yields 

V 4 P(x,y) = < V 4 G (i) (x,y) > [M 2 ] _1 {P. } 

i=l , 2 , . . . ,36 (2.18) 

where (V 4 G^ (x,y) ) [M 2 ] ^ is evaluated for any desired 

nodal point. The result is a molecule for the biharmonic 

operator with the desired specified boundary conditions. 

2 . Operators With Second Partial Derivatj e Boundary 
Conditions 

A biharmonic operator for a region with the harmonic 
operator specified along a part of the boundary is often 
required for the solution of a boundary value problem. The 
required molecule for the biharmonic operator with such speci- 
fied boundary conditions is developed using polynomial (2.2) 

and the grid of Figure (2) . The applicable equations for 
determining the constants of the polynomial are 

P(x,y) = <G^(x,y)> (D- } 

i=l , 2 , . . . , 36 (2.19) 

V 2 P (x,y) = < V 2 G (i) (x,y) > {D • } 

i=l,2,...,36 (2.20) 

from which the system of equations becomes 

{P,} = [M 3 ] {D ± } 

i=l , 2 , . . . , 36 (2.21) 

where 

T 

( p ± } = < p (Xj ,Yj ) | V 2 P (X k ,y k ) > 
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[M 3 ] = < G (i > (Xj ,yj ) > 



30 



/ • • • / 



< vV 1 * (x k ,y k )> 



k=31 , 32 



36 



i=l , 2 



/ • • • / 



36 



Solving for the unknown constants and operating on the 
resulting polynomial with the biharmonic operator yields the 
desired relationship for determining the biharmonic molecule. 
The relationship for the operator with the specified boundary 
condition is 

V 4 P(x,y) = < V 4 G ^ (x,y) > [M 3 ] 1 {D- } 

i=l,2,...,36 (2.22) 

For boundary conditions requiring the second partial 
with respect to y along one edge and the second partial with 
respect to x along a perpendicular edge, polynomial (2.2) and 
the grid of Figure (3) is used in the development of the 
applicable biharmonic molecule. The equations required to 
determine the appropriate constants for the interpolating 
polynomial are 



P(x,y) 



( G (l) (x,y) > {Dj.} 



i=l , 2 



36 



(2.23) 



P 



>yy 



(x,y) = < G (l) (x,y)> 




ryy 



/ • • • / 



36 



(2.24) 
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P ,xx ( x ' y) 



(2.25) 



V 



P 



P 



'P(x,y) 



,xxy 



(x,y) 



/Yyx 



(x,y) 



<G , ( xi (x ' Y)> {D i } 



i=l , 2 , . . . , 36 



< V 2 G (l) (x,y ) > {D.^ 



i=l, 2 , . . . , 36 



<G (l) (x # y)> {D.} 

' XX V 1 1=1,2,. ..,36 



<G(i) (x,y)> {D i } < 

'ix* 2 

J- -L / / • • • / 



36 



(2.26) 

(2.27) 

(2.28) 



Conditions expressed by equations (2.26), (2.27), and (2.28) 

are required to obtain sufficient independent equations with 
which to determine the thirty-six unknown constants. The 
values for these additional nodal boundary conditions are 
handled as additional unknowns. The resulting set of equations 
for determining the unknown constants of the required polynomial 
is 



{P^ = [M 4 ] {D ± } 



i=l , 2 , . . . , 36 



(2.29) 



where 



{P i> = < p < x j'Vj> ! p ,yy <x k' y k ) ! P ,xx (x £' y l > 



V 2p (x 26 /y 2 6), p ,xxy * x 35' y 35*> P ,y yx <x 36 ,y 36* ^ 
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[m 4 ] = 



< G^> (x,,y,) > 

J J j=l,2,..., 



25 



< G / ( w (x k^k ) > 

,yY k=29 , 28 , . . . , 30 



< > 



&=31 , 32 , . . . ,34 



(a) 



< v 2 g (1) (* 26 ,y 26 > > 



< G , < xiy (X 35' y 35 )> 



< G ,*yyx (x 36 '^36^ 

i=l , 2 , . . . ,36_ 
The final form of the desired finite difference 
biharmonic operator is 



V 4 P(x,y) = < V 4 G ^ (x,y) > [M 4 ] “ 1 {D i > 



i=l , 2 , . . . , 36 (2.30) 

3 . Operators With Mixed Boundary Conditions 

The development of a biharmonic molecule with boundary 
conditions of the first partial with respect to y along one edge 



( a ^It should be noted that points 26, 35, and 36 all have 
the same geometrical co-ordinates. 
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and the harmonic condition along a perpendicular edge utilizes 
interpolating polynomial (2.2) and the grid of Figure (3). The 
applicable equations for determining the required constants 
for the interpolating polynomial are 



P (x,y) 

P v (x,y) 
V 2 P(x,y) 

P ,xy (x '^ 



< G (i) (x,y) > 

< (x,y) > 
r y 

( V^G ^ (x, y ) 

< G ( j> (x,y) > 

r ^y 



{D i>. 

1=1 , 2 , . . . , 36 



{ Di } 

1=1 ,2 , ,36 



> (D.) 

i=l , 2 , . . . , 36 



{D.} 

1 i=l , 2 , ,36 



(2.31) 

(2.32) 



(2.33) 



(2.34) 



The resulting system of equations is 



{P i } = [M,] {D.} 

i=l ,2, . . . ,36 



(2.35) 



where 



<p.) = 



< P (x j,y j) 



IP 



,y (x k' 






I o 

I V Z P (x, 



.y t > 



IP 



, xy ^ x 36 



.y 36 )> 
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[m 5 ] = 



The final form of the desired finite difference 
biharmonic operator is 

V 4 P(x,y) = < V 4 G (l) (x,y) > [Me]" 1 {D.;} 

i=l , 2 , . . . , 36 (2.36) 

4 . Other Finite Difference Operators 

Certain physical boundary conditions include values 
of material properties in their evaluations. For example, in 
plate bending two such conditions are the equations for the 
moment and the shear along the free edge of a flat plate under 
some loading. The technique of combining the biharmonic 
operator and such boundary conditions into a single finite 
difference molecule is not practical in such cases. The reason 
being that in order to invert the matrices encountered in the 
development of the desired molecule, the material property 



< G (i) (x. ,y . ) > 

. 3 3 i=l,2,. ..,25 



< G ,y >( W> 



k=26 , 27 , . . . ,30 



< V 2 G (i) (x £ ,y £ ) > 



£=31,32, ... ,35 



< G Jxy (x 36' y 36 ^ ^ 



i=l , 2 , . . . ,36_ 
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would have to be specified before the inversion of matrix [M] . 
This would limit the application of the molecule to a specific 
material . 

To prevent this limitation, a molecule for the 
desired boundary condition alone is developed using the con- 
sistent interpolating polynomial (2.1). The technique is one 
of developing individual molecules for the partial derivatives 
involved in the boundary condition. The individual molecules 
are then combined with the applicable physical properties to 
form a general molecule which will represent the boundary 
condition . 

An example of this approach is illustrated in the 
development of molecules for the moment and the shear along 
the free edge of a flat plate. The governing equations are 
(Ref. 5) 

M n = D ( p ,nn< x 'Y) + V 1 p ,tt( x 'Y)) = 0 ( 2 . 31 ) 

and 

v = D(I >nnn (x 'Y) + (2-y) p ,ntt (x 'Y)) = 0 (2.38) 

where n is the normal co-ordinate, t the tangential co-ordinate, 
and y is Poisson's ratio and D the flexural rigidity of the plate. 

Molecules for each of the derivatives in equations 
(2.37) and (2.38) are developed by using the interpolating 
polynomial in the form given by (2.5). Operating on the 
interpolating polynomial with the required differential oper- 
ators yields the following relationships 
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nn (x ' y) 


= ( 


F ( ^(x,y) > 
, nn 


[M]" 1 


{P i } 

i=l / 2 , . 


. . ,25 


(2.39) 


tt (x,y) 


= ( 


F ( ^! (x,y) > 
/ tt 


[M] -1 


{P i } 

i=l, 2, . 


. . ,25 


(2.40) 


nnn (x ' y) 


= < 


F , < nn„ <x ' y)> 


[M]" 1 


{P i } . 

i = l / 2 , . 


. . ,25 


(2.41) 


ttn (x ' y) 


= ( 


p , < ttn< x 'y> > 


EM) - 1 


{p i } . , „ 







Evaluation of the above relationships at a nodal point along 
the boundary of the grid given in Figure (1) yields the 
required molecule for each derivative. 

The molecules are then combined to form the following 
finite difference approximations for equations (2.37) and 
(2.38) . 



l^L = D( ( (x,y) > [M]" 1 + y <F^)(x,y)> [M] 1 ) {P^} 

11 ,nn t u. l. j- 

(2.43) 

V = D( < F { i> (x,y)> [m] -1 + (2-y) <F ( j;|(x,y)> [M]" 1 {p.} 

z nnn / u i-n 

(2.44) 



D. ERROR OF APPROXIMATION 

The molecules developed in this work are based on approxi- 
mating a two dimensional surface wit’ a polynomial. Before the 
results can be effectively used, a bound on the error entailed 
in this approximation is essential. 
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The error of approximation for a finite difference 
biharmonic operator for some nodal point (m) is defined as 




(APPROX.) (2.45) 



The approximate value of the biharmonic operator at nodal point 
(m) is determined from the appropriate computational molecule 
and the nodal values associated with the molecule. 

An evaluation of the required nodal values appearing in 
the approximate operator, can be obtained from a Taylor Series 
Expansion in two variables. If a Taylor Series Expansion is 
centered on nodal point (m) , values for all other nodal points 
in the grid may be expressed in terms of the nodal value at 
(m) , its partial derivatives, and the grid spacing. Values for 
the nodal derivatives specified by the boundary conditions 
are also expressed in terms of the nodal point (m) by using 
the respective derivatives of the Taylor Series Expansion. 

Consequent evaluation of the required nodal values in 
terms of nodal point (m) and substitution into relationship 
(2.45) yields the following relationship for the error of 
approximation 



n . +4 



E t - Z Ei h 
i=l 



00 




m 




(n i +4-c i ) 



(2.46) 



n . =2 , 3 

l 



r * * * r 



oo 




00 



E i =CONSTANT 
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For the purpose of this work the above series was truncated 
after the first non-zero term in order to obtain a bound on 
the error for the respective approximate biharmonic operator. 
After truncation, the error for all of the approximate bihar- 
monic operators developed in this work was found to be propor- 
tional to at least the square of the mesh size (h) and a sixth 
partial derivative with respect to x and y. 

With the error proportional to at least the square of the 
mesh size (h) used for the finite difference approximation, 
resulting solutions should rapidly converge to the true solu- 
tion as the finite difference mesh is refined. The true 
solution should be obtained as the mesh size tends toward zero. 

To obtain an error proportional to at least the square 
of the mesh size spacing in all cases required the use of 
polynomial (2.2) for the development of all biharmonic molecules 
with boundary conditions included. 
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III. RESULTS 



The molecules for the various finite difference operators 
developed by the author in this work are presented in Tables 
(1) through (8) . 

Table (1) represents the molecule for the general bihar- 
monic operator evaluated for the center nodal point. This 
molecule is applicable for all nodal points in a mesh which are 
at least two nodal points removed from all boundaries. 

Table (2) represents the molecules for nodal points (1, 

3) and (2,3) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial with 
respect to y along an edge. Such a case would be the solution 
of a flat plate bending problem with a built-in edge restraint. 

Table (3) represents the molecules for nodal points (1,3), 
(2,3), and (1,2) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial along two 
mutual perpendicular edges. Such a case might be the solution 
for the deflection of a flat plate under uniform pressure with 
all edges built-in. 

Table (4) represents the molecules for nodal points (1,3) 
and (2,3) to be used for the biharmonic operator when boundary 
conditions specified by the harmonic operator are required. 

Such a boundary condition might be the moment for a flat plate 
with an edge simply supported. 

Table (5) represents the molecules for nodal points (1,3), 
(2,3), and (1,2) to be used to numerically evaluate the bihar- 
monic operator when boundary conditions specify the second 
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partial with, respect to x along one edge of a boundary and 
the second partial with respect to y along another edge. Such 
a condition would exist in solving a plate bending problem 
when all edges are simply supported. 

Table (6) represents the molecules for nodal points (1,3), 
(2,3), and (1,2) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial with respect 
to y along one boundary and the harmonic operator along an 
adjacent boundary. Such a case would be for a flat plate 
with one edge built-in and an adjacent edge simply supported. 

Tables (7) and (8) represent the molecules for nodal 
points (1,0) and (2,0) to be used to numerically evaluate 
the boundary conditions governed by equations (2.37) and 
(2.38) . 

Below each molecule a common factor is given by which the 
tabulated entries of the molecule must be multiplied. This 
common factor was generated by converting the original compu- 
tational values which were in terms of repeating decimal 
fractions to rational fractions. 

Each molecule is presented in a form similar to the grid 
(Figures 1, 2, and 3) from which it was developed. Consequently 
each constant in a molecule applies to the respective nodal 
deflection or boundary condition value in the grid for which 
it was developed. 

An estimate of the error of approximation is given below 
each molecule. 
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Table 7 - Finite Difference Molecule For 
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Table 8 - Finite Difference Molecule For 




IV. DISCUSSION 



A INTEGRITY OF DERIVED FINITE DIFFERENCE OPERATORS 

Confidence in the integrity of the finite difference 
operators developed in this work was established by comparing 
the solutions for the deflection of a flat plate under uniform 
pressure by the finite difference method with the exact solu- 
tion given by Timoshenko (Ref. 6) and Roark (Ref. 7). 

Each problem studied was solved using the molecules 
developed in this work and general finite difference techniques. 
Mesh sizes of L/4 , L/6 , L/8, L/10, L/12, where L is the width 
of the square plate studied, were used in successive solutions 
to the problem. The solution for the maximum deflection, 

^max f t * le P^ ate was use( ^ as a comparison. The value obtained 
from successive refinement of mesh spacing for each problem was 
plotted versus 1/N , where N is the number of divisions on one 
side of the plate. The resulting plots were then extrapolated 
to obtain the solution for a mesh having an infinite number of 
nodal points or the condition where the mesh spacing approaches 
zero. This extrapolated value is taken as the final solution 
for each problem and is compared with the exact solution. The 
above plots also indicate the convergence tendency of the 
solutions obtained with the molecules developed. 

The first problem solved was one of a square plate under 
uniform pressure with all edges built-in. The problem was 
solved using the molecules developed in this work, (Present 
Method) , and by means of the standard biharmonic molecule 
found in the literature and fictitious points, (Classical 



47 



Method) . The results obtained from these methods are pre- 
sented in Table (9). Figure C4) is the plot of successive 
solutions and illustrates the convergence trend for each 
method . 

A brief description is given in Appendix A of the solution 
process for the problems studied in this work. The solution 
process given is the one for the Present Method. 

The second problem solved was a square plate under uniform 
pressure with all edges simply supported. The results for the 
maximum deflection of the plate obtained from both methods of 
solution are presented in Table (10) . The plot of the results 
is given in Figure (5) . 

The third problem solved was a square plate under uniform 
pressure with two opposite edges built-in and the other two 
edges simply supported. A solution using the Classical Method, 
although possible, was not performed. The results for the 
maximum deflection at the center of the plate are given in 
Table (11) and the convergence plot by Figure (6). 

The last problem solved was a square plate under uniform 
pressure with two opposite edges simply supported, one edge 
built-in and the remaining edge free. The Classical Method 
is not feasible for the solution of this problem since the 
relationship for fictitious points outside the free edge is 
unknown. The use of a L/4 mesh spacing was also determined to 
be too coarse for a valid solution of the problem. The reason 
being that equations in addition to the governing equilibrium 
equation must be written to account for the boundary conditions 
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MESH 

SIZE 

SPACING 


NO. OF 
MESH 

DIVISIONS 


PRESENT 

METHOD 


CLASSICAL 

METHOD 


EXACT 
SOLN. 
REF. 6 


h 


N 


xio -3 pl 4 /d 


xio -3 pl 4 /d 


xio -3 pl 4 /d 


L/4 


4 


1.0455 


1.7995 




L/6 


6 


1.1293 


1.5343 




L/8 


8 


1.1813 


1.4245 




L/10 


10 


1.2095 


1.3696 




L/12 


12 


1.2258 


1.3388 




0 


00 


1.2628 


1.2688 


1.2637 


Table 9 - Maximum Deflection of a Square 


Plate Built- 




in On 


All Edges 


Under Uniform 


Pressure 


MESH 


NO. OF 


PRESENT 


CLASSICAL 


EXACT 


SIZE 

SPACING 


MESH 

DIVISIONS 


METHOD 


METHOD 


SOLN. 
REF. 6 


h 


N 


xlO -3 PL 4 /D 


xlO -3 PL 4 /D 


xlO -3 PL 4 /D 


L/4 


4 


3.6555 


4.0283 




L/6 


6 


3.8405 


4.0483 




L/8 


8 


3.9265 


4.0547 




L/10 


10 


3.9719 


4.0575 




L/12 


12 


3.9979 


4.0590 




0 


oo 


4.0570 


4.0630 


4.0567 



Table 10 - Maximum Deflection For Square Plate Supported 
On All Edges Under Uniform Pressure 
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xlO 3 PL 4 /D 




FIGURE 4 - Convergence Plot For Plate With 
All Edges Built-in 
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xlO 3 PL 4 /D 




FIGURE 5 - Convergence Plot For Plate With 
All Edges Simply Supported 
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MESH 

SIZE 

SPACING 


NO. OF 
MESII 

DIVISIONS 


PRESENT 

METHOD 


EXACT 
SOLN. 
REF. 7 


h 


N 


xlO -3 PL 4 /D 


xio -3 pl 4 /d 


L/4 


4 


1.6563 




L/6 


6 


1.7625 




L/8 


8 


1.8213 




L/10 


10 


1.8530 




0 


00 


1.9094 


1.9230 


Table 


11 - Maximum Deflection For Square Plate With 
Two Opposite Edges Built-in And The Other 
Two Simply Supported Under Uniform Pressure 


MESH 

SIZE 

SPACING 


NO. OF 
MESH 

DIVISIONS 


PRESENT 

METHOD 


EXACT 
SOLN. 
REF. 6 


h 


N 


-2 4 

XlO PL /D 


-2 4 

xlO PL /D 


L/6 


6 


0.96672 




L/8 


8 


1.0017 




L/10 


10 


1.0320 




L/12 


12 


1.0533 




0 


oo 


1.1017 


1.1263 



Table 12 - Maximum • Deflection For Square Plate With 

One Edge Built-in, One Edge Simply Supported, 
One Edge Free And One Edge Simply Supported 
Under Uniform Pressure 
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A PRESENT METHOD 



xlO 3 PL 4 /D 
1.94 r 




1.60 ' 1 * 1 1 1 1 1 1 

0.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 

1/N 2 

FIGURE 6 - Convergence Plot For Plate Built-in 
On Two Opposite Edges and Simply 
Supported On The Other Two Edges 
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xlO 2 PL 4 /D 




FIGURE 7 - Convergence Plot For Plat With 
Two Edges Simply Supported, One 
Edge Built-in And One Edge Free 
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along the free edge, since these conditions were not incorpor- 
ated into the finite difference biharmonic operator. The 
coarse grid given by a L/4 mesh does not provide sufficient 
nodal points to write equations to establish complete equil- 
ibrium over the entire plate. The requirement of additional 
equations for this particular problem is made clearer in the 
solution process in Appendix A. The results for the maximum 
deflection of the plate are given in Table (12), and the 
convergence plot is illustrated by Figure (7) . 

B. CONCLUSIONS 

Sets of mathematical molecules for the numerical approx- 
imation by finite differences of the biharmonic operator with 
various boundary conditions included were developed, as well 
as molecules for other differential operators. One consistent 
form of an interpolating polynomial was used to approximate 
the required two dimensional surface for the construction of 
all finite difference operators. 

Solutions for various plate bending problems obtained by 
using the developed finite difference operators are in 
excellent agreement with the exact solutions found in the 
literature . 

The convergence plots illustrate a straight convergence 
rate without oscillations for mesh spacings greater than L/6 
for all problems where the biharmonic operator and boundary 
conditions may be combined into one molecule. In Figures (4) 
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and (5 ) , solutions obtained by the present method are compared 
to the classical solutions with fictitious points. In one 
case the rate of convergence of the method is better than the 
classical method, in the other it is not. However the present 
method always gave a lower bound which is not always true for 
the classical method as shown in the two figures. 

The major advantage of the molecules developed is that 
they may be used to solve boundary values problems without the 
use of fictitious nodal points. Consequently, solutions to 
problems such as the plate with free edge boundary condition 
can now be solved using finite difference operators. 

One final conclusion concerning the overall aspect of 
the technique proposed should be made. Although the technique 
for developing finite difference operators proposed in this 
work has certain advantages over the classical method in the 
literature, it should be realised that other numerical methods 
exist for solving boundary value problems. The FINITE ELEMENT 
method has been applied very successfully and is covered in the 
literature very thoroughly, for example Refs. (8, 9 and 10). 

The method of finite differences using the consistent 
operators developed in this work is not superior to that of 
finite element methods for the solution of problems that may 
be formulated for a finite element scheme of solution, such as 
plate bending problems. However, difference operators may be 
developed, using the technique proposed in this work, to 
obtain an approximate solution to any problem for which the 
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governing equation is known. Thus they are applicable in the 
solution of problems that are not easily formulated into a 
finite element scheme of solution. 



s 
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APPENDIX A 

BRIEF DESCRIPTION OF SOLUTION PROCESS 

The solution of the plate bending problems by means of 
the finite difference operators developed in this work was 
performed by solving the governing equation V^w=p/D where 
w=w(x,y) is the displacement component perpendicular to the 
middle plane (x,y) of the plate. 

For the plate with all edges built-in, a mesh spacing 
of L/4 was used for the first solution. The applicable mesh 
is given below in Figure (8) . 




FIGURE 8 - L/4 Finite Difference Mesh For A 

Square Plate With All Edges Built- 
in Or Simply Supported 

Due to symmetry and the fact that w=0 along the boundaries , 
there are only three unknown values of deflections as shown 
in Figure (8) . For each unknown nodal deflection, the 
governing equation for plate equilibrium is written with the 
biharmonic operator being approximated by the appropriate 
difference operator developed in this work. The resulting 
set of simultaneous equations are then solved for the unknown 
deflections. For this problem the following molecules were 
used to express the governing equation 
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1 . 



Molecule from Table (3) for co-ordinate (1,3) 
applied at nodal point one. 

2. Molecule from Table (3) for co-ordinate (2,3) 
applied at nodal point two. 

3. Molecule from Table (1) applied at nodal point 
three . 

The derivatives required for the nodal points along the bound- 
aries are set to zero, thus establishing zero slope along the 
built-in edges. 

For the plate with all edges simply supported, the same 
mesh (Fig. (8)) is applicable for the first finite difference 
solution . 

The three unknown deflections are determined in this case 
by evaluating the governing equation with the following 
molecules 

1. Molecule from Table (5) for co-ordinate (1,3) 
applied at nodal point one. 

2. Molecule from Table (5) for co-ordinate (2,3) 
applied at nodal point two. 

3. Molecule from Table (1) applied at nodal point 
three . 

The derivatives required for the nodal points along the 
boundaries are set to zero, establishing zero moment for the 
simply supported restraints. 

For a plate with two opposite edges built-in and two 
edges simply supported, the L/4 mesh shown in Fig. (9) is 
applicable 
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FIGURE 9 - L/4 Finite Difference Mesh For A 
Square Plate With Top And Bottom 
Edges Built-in And Side Edges 
Simply Supported 

In this case, after applying symmetry and the fact of zero 
deflection along the edges, there are four unknowns to be 
determined. The required equations are written using the 
following molecules 

1. Molecule from Table (6) for co-ordinate (1,3) 

applied at nodal point one. 

2. Molecule from Table (6) for co-ordinate (2,3) 

applied at nodal point two. 

3. Molecule from Table (6) for co-ordinate (1,2) 

applied at nodal point three. 

4. Molecule from Table (1) applied at nodal point 

four . 

As before, the values of the derivatives specified by the 
boundary conditions are set to zero. 

For the plate with built-in, simply supported, free and 
simply supported edges, a mesh spacing of L/6 was used for 
the first finite difference solution as shown in Fig. (10). 
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13 
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13 




16 
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18 


17 


16 





FIGURE 10 - L/6 Finite Difference Mesh For A 
Square Plate With Top Edge Built- 
in, Sides Simply Supported, And 
Bottom Edge Free 

After applying symmetry and the fact that the deflection is 
zero on all boundaries except the free one, eighteen unknown 
deflections remain in the L/6 mesh. These unknowns are 
determined by solving the governing equation for nodal points 
1 through 12 and the boundary condition Equations of shear 
and moment along a free edge for nodal points 16, 17, and 18. 
This requirement of using the nodal points along the free 
edge to specify the applicable boundary conditions at this 
location reduces the number of nodal points for which the 
governing equilibrium equation may be written. Consequently, 
an L/4 grid spacing was found to be too coarse for a valid 
solution of the problem. 

The required eighteen equations are written using the 
following finite difference molecules; 

1. Molecule from Table (6) for co-ordinate (1, 

3) applied at nodal point one. 

2. Molecule from Table (6) for co-ordinate (2, 

3) applied at nodal point two. 

3. Molecule from Table (2) for co-ordinate (2, 

3) applied at nodal point three. 
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4. Molecule from Table (6) for co-ordinate (1, 

2) applied at nodal point four. 

5,6,8, Molecule from Table (1) applied at nodal 
9,11, points 5,6,8,9,11, and 12. 

12 . 

7,10. Molecule from Table (4) for co-ordinate (2, 

3) applied at nodal points seven and ten. 

In this case, the molecule is laid on the 
finite difference mesh in the y direction 
so that the specified boundary condition 
constructed in the operator is along the 
simply supported edge. 

13. Molecule from Table (7) for co-ordinate (1, 

0) is applied at nodal point sixteen. 

14,15. Molecule from Table (7) for co-ordinate (2, 

0) is applied at nodal points seventeen and 

eighteen . 

16. Molecule from Table (8) for co-ordinate (1, 

0) is applied at nodal point sixteen. 

17,18. Molecule from Table (8) for co-ordinate (2, 

0) is applied at nodal points seventeen and 

eighteen . 

Values for the derivatives along the built-in and simply 
supported edges required in the solution of this problem are 
set to zero to account for the specified boundary conditions. 
In order to obtain sufficient values for the deflection 



at a particular nodal point for the construction of a con- 
vergence plot, successive solutions are performed using a 
refinement in the mesh spacing for each solution. 
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